#!usr/bin/perl -w
use strict;

#under the directory of yzhou@biocluster:~/projects/fungal_TE/ful_len_finding/blast2seq
#run like this: perl ../TIR_finding/searchforTIR.pl -t 5_vs_4_ncbi.blastn.filterd.tab -db 5_Nc_10.fasta
#the output file is 5_Nc_10.full.fasta in the directory of Full

use Bio::DB::Fasta;
use Bio::SeqIO;
use Getopt::Long;
#use Bio::Perl;

my $db;
my $table;

GetOptions('d|db:s' => \$db,
			't|table:s' => \$table,
			);
			
			
open(TAB, "<$table")  || die "cannot open file $table: $!\n";
open(DATA, "<$db")  || die "cannot open file $table: $!\n";

my $tmp = $db;
$tmp =~ s/\.fasta//;
my $tmp1 = $tmp;
my $output = "./Full_TcMariner/$tmp.full.fasta";



my $qstart;
my $qend;
my $seqID; #we use query firstly

#record the alignment region of the query sequence from the blastn result table 
while(<TAB>)
{
	next if /^#/;
	my @row = split(/\s+/,$_);
	$qstart = $row[6];
	$qend = $row[7];
	$seqID = $row[0];
}
print "qstart is $qstart and qend is $qend\n";

my @leftend;
my @rightend;
my $rev_com_right;
my @EX_LF;
my @EX_RT;
my @full_seq;
my $full;
my @full_TSD;
my $full_TSD;
while(<DATA>)
{
	next if /^>/;
	my $count;
	my @row = split(//,$_);
	
	#this extracts the left flanking region of the sequence
	#print "the $qstart character is:";
	#print "$row[$qstart-1]\n";
	for(my $i=-10; $i < 100; $i++) # we have to check several nucleotides forward to make sure not miss TSD
	{
		push(@leftend, $row[$qstart-1+$i]);
	}
	print "the leftend seq is @leftend\n";
	
	#this extracts the right flanking region of the sequence
	for(my $i=-10; $i < 100; $i++) # this stores the reversed sequence of the rightend
	{
		push(@rightend, $row[$qend-1-$i]);
	}
	my $tmp = join("",@rightend) ;
	$rev_com_right = com($tmp); # this function gives the complement sequence
	@rightend = split(//, $rev_com_right);
	print "the reversed-complement rightend seq is @rightend\n";



#######################################matching section beginning###############################################	
	#this finds the longest common substring in the left and right end
	my $longest_com_sub;
	my $length=0;
	my $l = join("",@leftend);
	print "left string is $l\n";
	my $r = join("",@rightend);
	print "right string is $r\n";
	for(my $start=0; $start <= length($l)-3;$start++)
	{
		for(my $len = 3; $len <=length($l)-$start; $len++)
		{
			my $sub = substr($l,$start,$len);
			my $pos = 1+ index($r,$sub);
			if($pos>=1 && length($sub) >= $length) 
			{
				$longest_com_sub = $sub;
				$length = length($sub);
			} 
			#printf"str1: %d		str2: %d	'%s'\n",$start, $pos, $sub if $pos;
		}
	}

	printf"length: %d	'%s'\n",$length, $longest_com_sub;
#####################################matching section end#######################################################
	
	#make sure the two neuclotides before are "TA" in leftend and "AT" in rightend
	my $l_pos = index($l,$longest_com_sub);
	my $r_pos = index($r,$longest_com_sub);
	printf "left: %d 		right: %d\n", $l_pos, $r_pos;
	printf"leftTSD: %s%s 		rightTSD: %s%s\n",$leftend[$l_pos-2],$leftend[$l_pos-1], $rightend[$r_pos-2], $rightend[$r_pos-1]; 


	#judge if a full element was found
	my @l_TSD;
	my @r_TSD;
	for(my $i = $l_pos-2; $i <= $l_pos-1; $i++)
	{
		push(@l_TSD,$leftend[$i]);
	}
	
	for(my $i = $r_pos-2; $i <= $r_pos-1; $i++)
	{
		push(@r_TSD,$rightend[$i]);
	}
	
	
	my $lTSD=join("",@l_TSD);
	my $rTSD=join("",@r_TSD);
   if($length >= 9 && $lTSD eq "TA" && $rTSD eq "AT")
  {
	open(NEW, ">$output") || die "cannot open file $output: $!\n";
	print NEW ">";
	print NEW $tmp1;
	print NEW "\n";
	#retreive the full length element
	for(my $i=$qstart-11+$l_pos+2; $i <= $qend+9-$r_pos-2; $i++) # we have to check several nucleotides forward to make sure not miss TSD
	{
		push(@full_seq, $row[$i]);
	}
	$full = join("",@full_seq);
	my $start_index=$qstart-11+$l_pos+2; # 2 is to exclude TSD but we can keep TSD for double-checking
	my $end_index=$qend+9-$r_pos-2;
		
	printf "fullstarts: %d 		full ends: %d\n", $start_index, $end_index;
	
	for(my $i=$qstart-11+$l_pos; $i <= $qend+9-$r_pos; $i++) # we have to check several nucleotides forward to make sure not miss TSD
	{
		push(@full_TSD, $row[$i]);
	}
	$full_TSD = join("",@full_TSD);
	#for double-checking
	print "the full seq plus TSD is is:\n$full_TSD\n";

	print NEW $full;
  }
  else {
	print "L_TSD:$lTSD\n";
	print "R_TSD:$rTSD\n";
	print "full Tc1/Mariners not found in $db!\n";}
}



#this function returns the complemented sequence for a sequence 
sub com
{
	my $seq = $_[0];
	$seq  =~ tr/ACGTacgt/TGCAtgca/;
    return $seq;
}
			
			

